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ABSTRACT 

"s'  A  coapatational  algor itha  is  givan  for  obtaining  asyaptotieally 
officiant  astiaatas  of  tha  unknown  eoaplax  aaplitodas  and  frequencies 
in  a  supariaposad  exponential  aodal  for  signals.  It  is  shown  that  the 
variance  eovariance  aatrix  of  these  astiaatas  are  asyaptotieally , sane 
as  that  for  tha  aaxiaoa  likelihood  estiaatas  and  thus  attain  the 

Craaer-Rao  lower  bound.  ' 
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1 .  INTRODUCTION 

We  consider  the  superimposed  exponential  signals  model 
iu>  t  iu  t 

yt  *  aie  +  ...  +  ape  P+ft*  *  =  (1.1) 

where  i  =  V-T,  a.  is  a  complex  amplitude,  w.  is  the  frequency  of  the 
s-th  signal  and  {e^}  is  a  sequence  of  iid  complex  random  variables 
such  that 

E(ft)  =  0,  Eire  *t)a  »  E(isi  et)a  = 

co vire  €t,  im  e^)  *  0  (1.2) 

where  re  and  im  indicate  real  and  imaginary  parts  of  a  complex 
number . 

The  problems  of  interest  are  the  estimation  of  the  unknown 
frequencies  assumed  to  be  different,  and  the  estimation 

of  the  unknown  complex  amplitudes  a1,...,ap. 

When  has  a  complex  normal  distribution,  the  MIJS’s  (maximum 

likelihood  estimators)  of  the  unknown  parameters  are  the  same  as  the 
non-linear  LSE’s  (least  squares  estimators)  obtained  by  minimizing 

N  p  iw  t  , 

E  |yt-  £  a  e  *  |a  (1.3) 

t*i  x  s-i  8 

with  respect  to  a8  and  «s,  s  *  1 . p.  Unfortunately,  there  is  no 

closed  form  solution  to  this  problem.  In  this  paper,  we  develop  a 
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cooputer  algorithm  which  will  provide  estimates  which  are  asymptotically 
as  efficient  as  the  MLE’s. 

Let  us  write  the  unknown  parameters 


re  a  *  (re  a, , . . . ,  re  ol) 

~  1  O 


(1.4) 


im  a  »  (im  a, , . . .  ,im  a) 

~  l  p 


(1.5) 


«  *  («j. ••••"») 


(1.6) 


and  denote  the  Fisher  information  matrix  for  all  the  parameters 
(1.4M1.6)  by  Fjj.  Further  let 


Tn  »  diag(Nl/aI  ,  lfI/aI  ,  Hs/il  ) 


where  Ip  is  the  identity  matrix  of  order  p  and  define 
A  =  diag(at , . . . ,ap) .  Then  we  have 


(  21. 


lin  T^FJIC1 
*t . —  n  n  n 

inw 


-3 


21. 


-im  A) 
re  A 


v  »  • 
-im  A  re  A  -*A  i 


(1 ./) 


From  (1.7),  we  expect  that  w  and  a  »  (re  a,  im  a)'  the  MUS’s  of  uj  and 
a  have  the  limiting  distributions 


N*/a(£-w)  — »  Nn(0,6<ra(A*A)"1) 

a#  m  P 


(1.8) 


N1/a(o-a)  M,_(0,(raV) 

<V  ip 


(1.9) 
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where 


f  ylp-t~(//»i  A)3(A*A)_1  -  jire  A  im  A)(A*A)-1' 

-  -jkre  A  im  A)(A*A)-1  Tl*¥re  A)3(A*A)-1 

We  show  that  the  estimators  we  propose  have  the  same  limiting 
distributions  as  (1.8)  and  (1.9)  and  are  thus  asymptotically  efficient, 
i.e.,  attain  the  lower  bounds  for  their  asymptotic  covariances. 

The  best  known  methods  of  estimation  for  the  frequencies  u1 . wg, 

like  the  modified  FBLP  (forward  and  backward  linear  prediction)  of  Tufts 
and  Kumaresan  (1982)  and  EVLP  (equi variation  linear  prediction) 
discussed  in  Bai,  Krishnaiah  and  Zhao  (1986)  and  Bao  (1988)  have  certain 
deficiencies.  The  modified  FBLP  estimates  are  not  consistent,  although 
simulation  results  support  their  validity  in  small  samples  when  the  SNR 
(signal  to  noise  ratio)  is  relatively  low.  Hie  EVLP  provides  estimates 
of  the  frequencies  which  are  asymptotically  normal  and  have  a 

convergence  rate  of  0(N~1^3).  However,  this  is  still  not  the  best 
possible.  In  the  next  section,  we  describe  the  EVLP  method  and  show  how 
the  ELVP  estimates  could  be  refined  to  produce  fully  efficient  estimates 
of  the  frequencies  u  , . . . ,w  with  the  best  possible  convergence  rate 

A  3 

of  0(N"3/3). 


2.  THE  EVLP  METHOD 


Suppose  that  the  vector  b  *  (bQfb1 , . . . ,bp)  is  such  that 


+  b.z  +  ...  +  b  zp  =  b_  8  (z-e  *). 


o  l 


p  s=i 


(2.1) 
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Then  for  any  n  >  p+1 


S  bj  u  =  ^  by  £  l  (2.! 

feo^n_k  k*o^  n_k 

where  the  right  hand  side  of  (2.2)  is  a  function  of  error  only.  The 
coefficients  b.  are  estiaated  by  aininizing 


s  lJVn-klJ 

n = p+i  k=o 


(2.3) 


subject  to  the  conditions  bQ  >  0  and  |b|  *  1.  Such  a  method  of 

estimation  is  known  as  the  EVLP  method.  It  aay  be  noted  that  in  the  LP 
and  FHLP  methods,  the  expression  (2.3)  is  ainiaized  subject  to  the 
condition  that  bQ  *  1.  [Unfortunately,  the  restriction  bQ  *  1  is  not 

sufficient  to  ensure  the  consistency  of  the  estiaates  of  the  ratios  of 
the  b.  coefficients.] 


Now  write 


A 


7 


rs 


a  iww  . * 


and  construct  the  (p+l)x(p+l)  aatrix 


(2.4) 


f»(7Ti).  (2.5) 

It  is  easily  seen  that  the  EVLP  estiaate  6  of  b  is  the  unit  eigen 
vector  with  a  non-nega.ive  first  eleaent  providing  the  smallest  eigen 
value  of  f.  We  use  £  to  construct  the  polynomial  equation 
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b0  +  hiz  *  • 


A  M 

•  +  bpzp  =  °. 


obtain  solutions  in  the  fora 


(2.6) 


(2.7) 


and  take  u>, , . . .  .ST  as  estimates  of  or.  It  is  shown  in  Bai, 

1  P  1  ’  P 

Krishnaiah  and  Zhao  (1986)  that  w  is  a  consistent  estimate  of  w  with 

|  /  M 

a  convergence  rate  of  Op(N  ). 


3.  THE  MAIN  THEOREM 


Let  Ss  be  an  estimate  of  ws,  s  3  l,...,p  and  compute 


A 


w. 


(3.1) 


where 


N  „  -i£f  t  N  -iw  t 

C„  3  E  y.(t  -  |r)e  and  D  3  E  yfe 

t»i  4  2  n  t»i  ‘ 

Then  we  have  the  following  theorem. 

Theorem.  Suppose  that  satisfies  the  conditions  (1.2), 
«8,  s  3  1 . p,  ai a  distinct,  and 
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2s-w9  =  VM_1  5)’  6  6  (0’»]’  8  =  1 . P-  (3.2) 

Then 

(i)  w.-u.  *  0„(N_1_I<5)  if  <S  <  7  (3.3) 

8  8  p  4 

(ii)  N3/a(w-w)  — ►  M  (O^U^A)"1)  if  (3.4) 

/v  /v  D  4 

where  A  =  diag(a1 . ap) . 

In  the  next  section,  we  show  that  starting  with  the  EVLP  estimates 
<3  of  w  and  using  the  formula  (3.1)  repeatedly,  we  arrive  at  fully 
efficient  estimates  of  wg  having  the  limiting  distribution  (3.4). 

Proof.  We  have 


N  -iw  t  p  N  i  (<*/  — a?  )  t  N  -iw  t 
E  y. e  8  =  £  a  E  e  B  8  +  E  e.e  8 

t=i  1  m*i  "  t=l  t»i 


A  l 


a  J  (N)  +  RW • 

m  m 


(3.5) 


It  is  easy  to  see  that 


J_(N) 


'Op(l)  if  ■  *  s 

N  i(w  -w*)n 
Ijf+Kw.-w.)  E  t  e  88 


s  s 


t*i 


N  +  0  (N1"*)  if 
P 


(3.6) 
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where  uG  €  (wa,w),  and 

3  3  3 


R(N)  = 


N  -iu  t 
E  £fe  8 
t=l  1 

N  -iw  t  L-l  [-i(w_ -u  ) ] j  N  •  -iw  t 


=  E  e,e  9  +  E 
t=i  c  j=i 

0[N(w_-wJ]L  N 


3 _ 3 


TT 


s _ s 


L! 


s  l*tl 

t=i 


E  e+tJe 

t=i  ' 


(3.7) 


where  | 9 |  <1  and  L  6  >  1.  From  (3.7)  computing  the  order  of  the 
terms  on  the  right  hand  side,  we  have 


mi  l-1  o_or(1+*)j)  j+i/a 

R(N)  =  0  (N1/3)+  E  -2 - - - 0  (N  )+0  (1) 


Op(Nl/2). 


(3.8) 


The  expressions  (3.5)-(3.8)  imply  that 


N  -iu  t  e 

E  yte  8  =  a  N(l40n(N"®)).  (3.9) 

t*i  1  »  P 


Similarly,  one  can  prove 


N  M  — i£?  t  N  M  -iw.t  c 

J1yt(t"S)*  -  Eet(t-|)m  8  -  ios(fj(l+Op(N“5))(w8-w8). 

(3.10) 


By  (3.9)  and  (3.10)  w«  obtain 
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w  =  2  +  if  im 

3  S  N2 


w.+O  (N“^)(w  -«  )  +  if  /"(a"1  E  £.  (t  -  f)e  8 
s  p  88  N3  8  t=i  1  * 


a  N(l+0  (N“^) 
s  p 

N 


-iw.t 


). 


(3.11) 


Then  the  theorem  follows  from  (3.11)  using  the  following  fact 


.  /«  N  »  -iw  t 

{(&)  1/1  S  €t(t  -f)e  8  .  s  =  1,2 . p) 


t=i 


;tvw  2' 


CN(0,<r2Ip). 


(3.12) 


4.  RECURSIVE  ALGORITHM  FOR  ESTIMATION 


Ve  start  with  a  consistent  estimate  of  w  and  improve  upon  it 

S 

step  by  step  by  a  recursive  algorithm.  The  m-th  stage  estimate  uim)  is 

3 

computed  from  the  (m-l)th  stage  estimate  by  the  formula 


£<m-i, 

8 


m  a  1,2,... 


(4.1) 


where 


"m 


D.  *  .!,V  ■ 

t*l 


N_ 


O.  -  -  f-u 

•■l 


N_  -ifc*-l»t 


Bx  S 


(4.2) 


(4.3) 
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Ve  apply  the  formula  (4.1)  repeatedly  choosing  N  suitably  at  each 

in 

stage. 

Step  1  with  m  =  1.  Choose  N  =  [N0,4]  and  w‘0)  =  2;  the  EVLP 

1  8  8 

estimate.  Note  that 

-  Op(K-,/J)  =  Op[M“I_1/4].  (4.4) 

Then  substituting  N  =  [ N° * 4 ]  and  w*0)  =  £5  in  the  formula  (4.1)  we 

l  8  8 

find  by  applying  the  result  (3.3)  of  the  main  theorem 

=  On(N71-I/a)  =  0  (N~°‘8) .  (4.5) 

8  8  p  1  p 

Step  2  with  m  =  2.  Choose  Na  *  [N°'4a]  and  using  Wg1>-Ws  *  Op(N-0’6) 

=  Op(N~1-1^4),  compute  by  the  formula  (4.1).  Again  by  the  result 

(3.3)  of  the  main  theorem 

w‘2,-«s  =  Op(N"1“l/2)  *  Op(N”°-7J).  (4.6) 

Steps  3  to  7.  Choosing  N3,...,N7  as  given  below  and  applying  the  main 
theorem  in  the  same  way,  we  have 
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N3  =  CN°  * 57] 

yielding 

A(  3 ) 

^S 

=  Op(N-0,84) 

N4  =  [N0-67] 

n 

A(  4  ) 

^8 

=  0p(N-‘) 

„  _  fNo . SO-j 

5 

n 

A  (  S  ) 

W1  ^S 

=  OpOT1-30) 

N8  =  [N°*9®] 

ft 

A(  6 ) 

u 

•o° 

* 

(4.7) 

Finally  we  take  N7  =  N 

A  »  y  » 

and  compute 

Now  applying  the 

result  (3.4) 

of  the  main  theorem,  we  have 


N*/3(w(7,-«j)  — >  M(0,6<r,(A*A)"1)  (4.8) 

which  is  the  same  as  that  of  the  MLS  given  in  (1.8)  and  shows  that  w(7) 
is  a  fully  asymptotically  efficient  estimate  of  u. 

Using  the  estimate  u/(7>  of  w,  a  is  estimated  by 

a  9  (4.9) 

where  A  is  a  nxp  matrix  whose  (r,s)-th  element  is  exp(irw‘7)) .  It  is 

8 

A 

easily  shown  that  a  has  the  same  limiting  distribution  as  that  of  the 
MLS  given  in  (1.9). 

Remark  1.  Note  that  =  N(I+0„ (N-1))  so  that  (4.9)  can  be  more 


8 imply  approximated  by 
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A 

a 


(4.10) 


Remark  2.  The  recursive  formula  (4.1)  for  N_  can  be  improved  by 
using  the  alternative  formula 


N-N. 


O'”' 


O'"-11  ♦ 

s 


m  C. 


i2 _  E  “  iM  rai 

"Xvn  tao  D;t 


(4.11) 


where 


N  a 

°.t  *  £  V.'"  *  '  ‘ 


N  -A 

m  - iw  ,n 

D.t  ■  n?/n*t«  ‘  ■ 


5.  SIMULATION  RESULTS 

In  order  to  examine  the  behaviour  of  the  statistics 

T,  =  N3/*(w*7,-w  ) 

1  8  5 


T,  »  Nl^3(re  a -re  a) 

2  5  8 

T,  =  N1/a(/m  o.  -  im  am ) 

5  5  5 

for  s  a  1,2,3,  as  N  increases,  the  following  simulation  study  was 
carried  out. 
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Using  the  model 


i«  t  iw  t  iw.t 

Yt  3  V  +  V  +  “3®  +  ft 

t  »  1 . N 

for  3  signals  with  the  true  values 


u.  «  2.1,  w,  a  2.9 

2  3 

im  a.  =  0,  is  1,2,3 

a1  *  5.0 


Wj  a  1.5, 
re  flj  »  5.2, 


independent  samples  of  sizes  varying  from  50  to  2500  at  intervals  of  50 
were  drawn  and  the  estimates  of  all  the  parameters  were  computed  as 
described  in  Section  4. 

Simulations  were  done  using  the  two  methods  given  'n  (4.1)  and 
(4.11).  The  results  show: 

(i)  When  the  sample  size  n  is  greater  than  300,  there  is  not  much 
difference  between  the  two  methods.  In  this  case  (4.1)  provides  a  simple 
method. 

(ii)  When  the  sample  size  is  less  than  300,  the  results  by  (4.1) 
appear  to  be  less  stable  than  those  by  (4.11).  Even  when  n  is  as  small 
as  50,  the  results  by  (4.11)  show  considerable  improvement  over  the  EVLP 
method. 

The  graphs  of  the  statistics  T1#  Ta  and  T3,  obtained  by  the  method 

(4. lx),  against  the  sample  size  N  are  shown  in  Figures  1,  2,  and  3 
respectively.  It  is  seen  from  the  graphs  that  stability  is  reached  at  a 
sample  size  of  the  order  of  250. 
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Flgure  3.  Graph  of  agalnat  sample  slza  N 


The  estiaatei  of  obtained  after  7  steps  of  refinement 

over  the  EVLP  estimates  are  given  in  the  following  Table.  The  estimates 
are  close  to  the  true  values  even  for  small  sample  sizes. 


TABLE.  Eatiaates  of  u/j,  «Jt  w,  obtained 

after  7  steps  of  refineaent  over 
the  EVLP  estiaates. 


N 

A 

wi 

A 

wa 

A 

50 

1.4902339 

2.0988951 

2.8970912 

150 

1.5000117 

2.0995676 

2.9003642 

250 

1.4997012 

2.1000414 

2.9000647 

350 

1.5000774 

2.1000476 

2.8999323 

450 

1.4999175 

2.1000282 

2.8999789 

550 

1.5001009 

2.1000981 

2.9001332 

650 

1.4999101 

2.0999894 

2.9001015 

750 

1.4999827 

2.1001022 

2.9000483 

850 

1.5000169 

2.1000005 

2.9000164 

950 

1.4999677 

2.0999751 

2.8999948 

1050 

1.4999554 

2.0999370 

2.9000816 

1150 

1.4999777 

2.1000294 

2.8999919 

1250 

1.4999834 

2.0999702 

2.9000048 

1350 

1.4999863 

2.1000342 

2.9000128 

1450 

1.4999595 

2.1000282 

2.9000266 

1550 

1.5000098 

2.1000278 

2.9000022 

1650 

1.4999914 

2.0999927 

2.8999936 

1750 

1.4999845 

2.0999933 

2.9000010 

1850 

1.5000072 

2.1000128 

2.9000106 

1950 

1.5000064 

2.0999905 

2.9000050 

2050 

1.5000024 

2.1000091 

2.8999912 

2150 

1.4999927 

2.0999981 

2.9000172 

2250 

1.5000053 

2.1000023 

2.8999934 

2350 

1.5000038 

2.0999954 

2.8999866 

2450 

1.5000119 

2.1000125 

2.9000161 
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signals. 
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A  computational  algorithm  is  given  for  obtaining  asymptotically  efficient 
estimates  of  the  unknown  complex  amplitudes  and  frequencies  in  a  superimposed 
exponential  model  for  signals.  It  is  shown  that  the  variance  covariance 
matrix  of  these  estimates  are  asymptotically  same  as  that  for  the  maximum 

likelihood  estimates  and  thus  attain  the  Craa^r-Rao  lower  bound. 
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